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Abstract 


The Refined Zigzag Theory is applied to the modeling of 
delaminations in laminated composites. The commonly used cohesive 
zone approach is adapted for use within a continuum mechanics model, 
and then used to predict the onset and propagation of delamination in 
five cross-ply composite beams. The resin-rich area between individual 
composite plies is modeled explicitly using thin, discrete layers with 
isotropic material properties. A damage model is applied to these resin- 
rich layers to enable tracking of delamination propagation. The 
displacement jump across the damaged interfacial resin layer is 
captured using the zigzag function of the Refined Zigzag Theory. The 
overall model predicts the initiation of delamination to within 8% 
compared to experimental results and the load drop after propagation is 
represented accurately. 
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1. Introduction 


Composite materials, owing to their favorable specific strength and specific stiffness properties, are a key 
technology for reducing the weight and improving performance of future aerospace vehicles. An 
important aspect in designing advanced composite structures is understanding how these materials fail 
under service loads. The failure behavior of laminated composite structures is often governed by complex 
interactions of multiple intralaminar failure mechanisms, e.g. matrix cracking and fiber fracture, as well 
as interlaminar damage such as delamination. The damage mechanism that influences the load 
redistribution from failed to unfailed material is typically direction-dependent and is often governed by 
geometric and/or constitutive non-linearities. Furthermore, advanced composites maintain considerable 
load-carrying capability past the initiation of the first failure site. This behavior arises because the 
failed/damaged material has the capability to unload and redistribute the internal loads to the remaining 
undamaged material, and because adjacent plies can create a constraining effect that localizes the damage. 
Under these circumstances, initial first-ply failure criteria such as the maximum strain, maximum stress, 
Tsai-Wu, Tsai-Hill or Flashin-Rotem criteria [1] may only provide conservative estimates of laminate 
strength. As a result, progressive damage models that model the degradation of material and the growth of 
damage until final failure have received much attention over the past 30 years [2]. 

For the sake of simplicity, the study of damage in composites is often split into separate analyses of 
intralaminar damage and delaminations, even though these mechanisms may interact. Material 
degradation models such as ply-discounting and continuum damage models are two methods for 
modeling intraply damage. Ply-discounting techniques use heuristic knock-down factors to degrade 
certain material properties of a ply when failure is detected using an interactive or non-interactive 
initiation criterion. Both the magnitude of the knock-down factor and the degraded engineering constants 
may be failure-mode dependent. In continuum-damage models, the material properties are gradually 
degraded as a function of an evolving internal state variable, such as strain, that models the accumulation 
of damage. Therefore, these techniques are less heuristic than ply-discounting techniques and may be 
derived using thermodynamically consistent relations [3]. 

Traditionally, the study of delaminations is split into two distinct phenomena: the initiation of the crack 
and the ensuing propagation, where the latter is based on formulae of a pre-existing crack. The virtual 
crack closure technique (VCCT) proposed by Rybicki and Kanninen [4] is often used for predicting 
delamination growth. However, in the finite element method (FEM), VCCT results are dependent on the 
mesh in front of the crack tip and a pre-defined location of a crack is required. Unfortunately, the exact 
location and size of this initial crack may be difficult to determine a priori. 

One approach to overcome the limitations of VCCT is the use of interfacial decohesion elements. The aim 
of the decohesion elements is to model the quasi-brittle process zone that appears in the wake of the crack 
tip. Usually, propagation of a crack does not cause immediate separation of two layers, whereby the 
residual cohesion may occur due to fiber bridging across the crack, or sliding friction between the 
surfaces. A cohesive law relates the three tractions acting on the two fracture surfaces with the associated 
interfacial displacement jumps [5]. For example, a cohesive law for Mode 1 relates the direct normal 
traction to normal displacement jumps, while the laws for Mode II and Mode III relate the shear tractions 
to in-plane displacement jumps. In the FEM, the cohesive zone is most often implemented using user- 
defined cohesive elements to represent the inter-ply resin that separates the actual plies modeled using full 
three-dimensional (3-D) elements. However, the high computational cost of non-linear 3-D finite element 
models restricts the use of cohesive elements to pre-determined failure sites and prohibits the use of these 
models for preliminary design studies. 

Higher-order equivalent single-layer models are an attractive compromise between reduced computational 
cost and accurate 3-D stress fields. In this respect, the Refined Zigzag Theory (RZT) developed by 
Tessler et al. [6-8] is a potential candidate as it has been shown to accurately model the stretching and 
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bending of laminated beams and plates, even when the material properties of individual plies vary by 
orders of magnitude [9], This feature makes the RZT applicable to explicitly modeling finite interfacial 
resin layers in pristine and degraded (damaged) states. Currently, only one attempt has been made to use 
RZT to model delamination damage in laminated composites by applying a continuum damage model to 
an interfacial resin layer [10]. In the current paper, the concept of cohesive zones is applied within the 
framework of RZT to model delaminations in beams. 


2. Approach 

In RZT, the displacement field of a beam is given in terms of four functional unknowns, 


u k (x,z ) = u 0 (x) + Z0(x) + (j) k (z)cp(x) 

( 1 ) 

u z (x) = w(x) 

( 2 ) 


where u 0 is a uniform axial (stretching) displacement, 6 is the bending rotation, (p is the zigzag rotation, 

and w is the transverse displacement. The zigzag function (j) k in ply k is derived from the ratio of the 

local transverse shear modulus G xz to the equivalent transverse shear modulus of the whole laminate (see 

[6]). The zigzag function models the occurrence of slope changes in the axial displacement field u k at 

layer interfaces that occur because of differences in the transverse shear moduli of adjacent layers 
(Figures 1-5). 

The displacement field assumption is used alongside the principle of virtual displacements to derive a set 
of variationally consistent governing field equations and boundary conditions [6] that are solved for the 
four functional unknowns. In RZT, the axial stress field <7 X is derived from the kinematic equations and 

constitutive relations. The transverse shear stress T xz and transverse normal stress <7. are derived in a 
post-processing step by integrating the axial stress in Cauchy equilibrium equations for greater accuracy. 

The present approach of modeling delaminations using RZT is best explained using an example. Consider 
a [O12/9O/O12] laminated beam, simply supported at both ends and loaded at mid-span by a transverse point 
load. As the Young's and shear moduli E x and G xz of the 90-degree ply are degraded, the displacement 

field transitions from purely linear to a pronounced zigzag shape through the thickness. When the 90- 
degree ply is fully damaged, two observations can be made regarding the laminate behavior: 

• The axial displacement is fully reversed throughout the O 12 sublaminates, i.e. tension and 
compression occur within each sublaminate. This means the sublaminates are bending 
independently as the 90-degree ply provides little structural connection (Figures 6-7). 

• Due to the pronounced shear deformation throughout the 90-degree ply, the O 12 sublaminates move 
relative to each other, i.e. there is a displacement jump between the bottom interface of the upper 
sublaminate and the top interface of the lower sublaminate (Figure 8). 

Thus, the approach presented herein is to insert layers of isotropic resin material between the composite 
plies and use a cohesive zone approach to model the interfacial degradation. In this manner, Mode II 
shearing fracture is modeled for beams, and both shearing Modes II and III are modeled for plates. Note 
that because RZT is based on an equivalent single-layer formulation, the material properties of the 
interfacial resin layers cannot be degraded absolutely to zero in order to maintain numerical stability 
when solving the governing equations. Therefore, in this paper, a fully damaged resin layer has material 
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properties that have been degraded by 12 orders of magnitude, resulting in negligible properties compared 
to the undamaged composite plies (Figures 6b and 7b). 

It is not possible to implement the cohesive zone concept directly within RZT because the commonly 
used cohesive element within FEM is not a full 3-D element based on continuum mechanics. Rather, it 
represents the interfacial surfaces that are created as the crack propagates and therefore defines the 
constitutive relation between the three surface tractions (one normal and two shearing) and the three 
components of the interfacial displacement jump (Figures 9-10). As there is no unique relation between 
stresses and displacements in 3-D elasticity, an arbitrary constant of proportionality K is defined for 

conventional cohesive elements used in 3-D FEM discretizations. Conversely, each ply within RZT 
requires a constitutive relation between the stress and strain tensors by using the elasticity tensor as a 
constant of proportionality. Thus, a pre-defined relation exists between stresses and strains based on the 
isotropic material properties of the resin. To remedy the inconsistency between these two approaches, the 
method proposed by Fan et al. [11] is adapted to RZT in this paper (Figure 11). The interfacial shearing 
displacement jump 5 k is derived from 



(3) 


u . — K 

where t is the thickness of the resin layer k and y xz is the average shear strain within the resin layer. 
This approach is equivalent to defining the initial penalty stiffness K of a cohesive element within a 3- 
D FEM discretization to be 


K 


p 



(4) 


as used by several authors in the literature [5], and where t is the thickness of the interfacial resin layer. 

Thus, in the framework presented herein, thin resin layers of 10pm thickness are modeled as explicit 
isotropic layers between composite plies, and the axial stress <j x , transverse shear stress r xz and 

transverse normal stress cr_ within this layer are calculated from the RZT governing equations. Damage 
initiation within a resin layer is defined when the quadratic failure criterion 


e = 





(5) 


exceeds unity, where T and S are the tensile and shear strengths of the resin, respectively. Macaulay 
brackets (. . .'j are used as compressive normal tractions do not contribute to the onset of delamination. At 

the point of failure initiation, the critical interfacial shearing displacement jump 5 k is calculated from the 
critical transverse shear stress T x , c within the resin layer, 



( 6 ) 


Thus, the damage variable d k at a point within a damaged resin layer is found using the usual bi-linear 
cohesive law [5], 
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( 7 ) 


.* a( s l -s‘) 
sLSe'-S') 

where S' mx is the maximum interfacial shearing displacement jump that has occurred during the load 
history and S k is the final interfacial displacement jump at which complete decohesion occurs. The value 
of S k is found from the fact that the area under the cohesive traction-displacement curve is equal to the 
critical surface fracture energy G Uc . Thus, 
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k 
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Once the damage variable d k at a specific mesh point within a resin layer is determined, all engineering 
constants within the elasticity tensor are degraded 

a k =(\-d k ]E k s k (9) 

As a result, the internal loads are redistributed within the laminate via two mechanisms. First, the 
reduction of the axial and transverse shear moduli causes a load redistribution to the surrounding material. 
Second, degradation of the transverse shear modulus G xz locally increases the zigzag function such that 

the sublaminates above and below a damaged interface can move relative to each. Without the added 
functionality of the zigzag function within RZT this mechanism would not occur. 

Accurate derivations of the transverse stresses T x _ and cr_ using Cauchy equilibrium equations in a post- 
processing step depend on accurate differentiation of the axial stress field <J X . In the classic weak-form 

FEM, the recovery of these gradients from the element shape functions can lead to the propagation of 
numerical noise. Discrepancies in the transverse stress profiles lead to errors in the damage initiation 
location. This point is important as small perturbations in the initial damage site can lead to different post- 
damage behavior. 

In this paper, the strong-form FEM based on elemental shape functions derived from the general 
differential quadrature method (GDQM) is used to discretize and solve the problem numerically. In 
GDQM, the governing differential equations are converted into algebraic ones by using a matrix 
differential operator that is derived from Lagrange polynomials (Figure 12). In this manner a derivative of 
a function can be replaced by a linear weighted sum of all functional values within the discretization 
domain, i.e. 


T/frl = Af'f(x j) / = 1,2,..., AT (10) 

OX 

where A‘" > is the n th derivative weighting matrix. Thus, the domain is discretized into finite GDQM 
elements that are based on higher-order polynomial shape functions giving an hp-type FEM. As the 
governing equations are solved in the strong form, derivatives of the axial stress cr x can be calculated 
accurately and efficiently using the GDQM weighting matrixes [12], and numerical noise in the 
derivation of transverse stresses r v _ and cr z is minimized. 
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Due to the degradation of material properties, a non-linear solver is required. In this work, an iterative- 
incremental Newton-Raphson solver was implemented in MATLAB® to solve the strong-form GDQM 
equations. Thus, the incremental equilibrium equation 

t+At,i ) ~~ f(t+M,i) ~ l+At,i ) (1 1 ) 

is solved until the residual vector R is lower than the convergence tolerance of 10 4 . In the equation 
above, ^(c{ +AM) ) is the tangent stiffness matrix at load step (t + At,i) in the i th iteration, and C k 
equals to the material stiffness matrices of the different plies in the laminate. The difference between the 
known external force vector F and internal force vector f gives the residual R , while A U( t+Al ( | is the 

incremental displacement vector at the i th iteration. Before damage occurs, the material stiffness matrices 
of all resin layers are equal to their pristine values C<J . However, upon damage initiation within a resin 

layer, the engineering constants that define C k are degraded at the damaged locations according to Eq. 9. 


3. Discussion and Conclusions 

The approach presented herein was used to predict the failure behavior of several cross-ply beams. For 
this purpose, the experimental results of the three-point bend test in the experimental study by Yuen et al. 
[13] were used for comparison purposes. In these tests, a single 90-degree ply was inserted at different 
locations throughout the thickness of a [O 24 ] laminate (Figure 13). The more compliant 90-degree ply acts 
as a weak interface to initiate delamination without the need for manufacturing a specimen with a 
preexisting crack. This load case is dominated by Mode II fracture and therefore presents a good 
validation study for the presented Mode II RZT damage model. 

In the present RZT model, interfacial damage is modeled by inserting thin 10pm isotropic resin layers on 
either side of the 90-degree ply using the resin properties given by Yuen et al. [13], The results in Figures 
14 and 15 indicate that the pre-damage behavior and post-damage load drop is modeled accurately by the 
RZT damage model. The load at which the first load drop occurs is predicted to lie within 8% of the 
experimental results for all cases considered (see the tables on p. 18). The location of the initial 
delamination site is predicted most accurately when the load pin used in the experimental studies is small 
(case S-40-1/2). This response happens because the contact area of the load pin was ignored in the 
computational model by applying a concentrated load. However, in all five cases, the size of the 
delamination after the first load drop is predicted accurately by the RZT model (Figures 16 and 17). 

In conclusion, the general trends of delamination initiation and propagation are modeled to adequate 
accuracy (~10% error) given the low computational cost of the underlying RZT equivalent single-layer 
beam model. The most accurate predictions, compared to the experimental results, occur when the load 
pin contact area in the experiments was small, as this contact area is not modeled explicitly herein. 

Currently, the resin layer is only modeled explicitly on either side of the 90° layer, as this is known to be 
the weak location. Further work needs to focus on placing resin layers within the 0° blocks to make sure 
that the model does not predict delamination here. 

The bending load case considered herein results in fracture behavior that is Mode II dominated. Thus, the 
present investigation can serve as a foundation for applying the methodology to end-notched flexure 
specimens and general laminates with off-axis fiber angles, and these load cases will be the topic of future 
work. The present approach should also be extended to the analysis of plates to analyze both Mode II and 
Mode III fracture. On the other hand, in-plane mode I-dominated load cases will be modeled less 
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accurately with the present framework, but this can be remedied by incoiporating transverse normal ZZ 
effects in the underlying RZT formulation. 

Finally, the cohesive zone approach presented herein is only one possible application of RZT to damage 
modeling. A large number of continuum damage models for interply and intraply damage have been 
published in the literature, and any one of these can be explored within the general framework presented 
herein. Perhaps, intraply and interply damage models can be combined within the presented framework to 
account for both modes of failure simultaneously. 
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Fig. 1 Arbitrary stacking sequence of A layers with notation of layer thicknesses and interface locations. The dashed 
line represents the continuous through-thickness deformation of CLA, whereas the solid, piecewise continuous line 

shows the ZZ displacement field. 



Fig. 2. The ZZ effect is driven by differences in transverse shear moduli of adjacent layers, e.g. 0° and 90° layers. As 
transverse shear stresses are continuous at layer interfaces, the transverse shear strains are different at the interfaces of 


distinct layers, and therefore cause a piecewise change in 
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Fig. 3 In the RZT beam theory, the Timoshe nk o beam displacement field is enhanced via a ZZ variable multiplied by a 
piecewise continuous ZZ function (j) k (z) . This function is derived by assuming that the transverse shear behavior of 
the laminate can be modeled as a system of springs-in-series with an equivalent spring stiffness G. 




Fig. 4 Differences between the RZT ZZ function in (a) and Murakami’s ZZ function in (b). Murakami’s ZZ function 
takes arbitrary values of ±1 at layer interfaces. The RZT ZZ function is zero at the top and bottom surfaces but is based 

on the actual material properties (adapted from [14]). 
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Fig. 5 The RZT ZZ function can model laminates where different layer have material properties that vary by orders of 
magnitude. For example, a sandwich beam with stiff face layers and heavily degraded central core, shows a full stress 
reversal in the outer layers when loaded in bending, i.e. the outer layers are both in compression and tension. CLA 

cannot account for this behavior. 
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(a) (b) 

Fig. 6. A [O12/9O/O12] laminated beam in three-point bending where the stiffness tensor of the central 90° layer is 
degraded using the damage parameter d. When the central layer is undamaged in (a), the displacement field is linear as 
predicted by CLA. Flowever, when the central layer is damaged in (b), a displacement jump occurs, which is modeled 

via the RZT ZZ function. 




(a) (b) 

Fig. 7 Similarly, the axial stress field throughout the laminate thickness changes from a classic solution when the 
central layer is undamaged in (a) to a non-classical stress field with stress reversal in the blocks of 0° when the central 

layer is damaged in (b). 
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Fig. 8 Schematic showing how the RZT beam model can be used to explicitly model the thin interfacial resin layers. 
The resin layer is modeled using isotropic material properties, and when undamaged does not influence the linear 
bending displacement field. When damaged, the resin layer allows the two sublaminates to slide (shear) over each 

other, and this is modeled using the ZZ function. 



Fig. 9 Schematic showing the concept of a brittle Mode I cohesive zone in a double cantilever beam test. The cohesive 
zone is modeled via a bi-linear traction-displacement law of fictitious initial stiffness. When a critical stress limit is 
reached the layer is damaged, and stress if offloaded onto other parts of the structure (adapted from [15]). 
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Deformed 3D element 



Deformed 3D element 


Fig. 10 Schematic diagram of a 3D cohesive finite element embedded between two 3D continuum finite elements. The 
cohesive element is not governed by the classic stress-strain constitutive tensor but by a pre-defined traction- 

displacement law (adapted from [16]). 



Mode II only for 



(a) 


(b) 


Fig. 1 1 Schematic showing the implementation of the cohesive zone damage model of Fan et al. [11]. Flere, the resin 
rich interfacial layer is not modeled using a classical cohesive element with a traction-displacement constitutive model, 
but rather using the classical stiffness tensor, i.e. a stress-strain law. Thus, the stiffness tensor is degraded directly when 
damage occurs. The displacement jump across the interface in (b) is derived from the thickness of the resin layer and 
the transverse shear strain across it. The continuum damage model approach means that this cohesive zone model can 
readily be implemented in a RZT beam model to analyze Mode II fracture (adapted from [11]). 


15 




NM 



Loading pin 



Fig. 13 The loading configuration of a laminated beam with two blocks of [0]i2 and one 90-layer modeled and 
experimentally tested by Yuen et al. [13]. The 90-layer is either placed in a neutral central position creating a 
[O12/9O/O12] laminate, or on the tension side for a [06/90 /Ois] laminate, or on the compression side for a [Ois/ 90/06] 
laminate. These configurations are denoted by V 2 beam, % beam and 3 /4 beam, respectively. Different length beams of 
40 mm, 48 mm and 56 mm were tested by Yuen et al. [13], and these results are treated as benchmarks herein. 
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Yuen et al. [13] 

Experimental Computational 



Fig. 14 Comparison of experimental results of 48 mm long beams by Yuen et al. [13] in (a) and the corresponding RZT 

beam computational results in (b). 
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Fig. 15 Comparison of experimental results of 40 mm and 56 mm long beams by Yuen et al. [13] in (a) and the 

corresponding RZT beam computational results in (b). 
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Specimen 

Crack initiation load (N) 
Experimental [13] Computational 

Error 

(O/o) 

48-1/4 

6.49 (0.37) 

6.60 

1.7% 

48-1/2 

6.22 (0.26) 

6.72 

8.0% 

48-3/4 

8.20 (0.09) 

8.60 

4.9% 


Loading pin contact area disregarded 
in model, i.e. a source of error in the 
predicted initiation location 


Fig. 1 6 Comparison of the initial delamination initiation site and final propagation length of experimental results by 
Yuen et al. [13] and the present RZT results. The table shows the percentage error in predicting the crack initiation 

load. 
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Crack initiation load (N) Error 

Experimental [13] Computational ( 0/, °) 

6.46 (0.26) 6.72 4 . 0 % 

6.59 (0.3) 6.72 2 . 0 % 

Fig. 17 Comparison of the initial delamination initiation site and final propagation length of experimental results by 
Yuen et al. [13] and the present RZT results. The table shows the percentage error in predicting the crack initiation 

load. 
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